# 02_NYMBY_AMCE.R
# Created: 2021-05-14 Taka-aki Asano
# Last Modified: 2024-02-16

# package
require("readr")
require("dplyr")
require("tidyr")
require("ggplot2")
require("ggpubr")
require("cjoint")
require("cregg")
theme_set(theme_minimal(base_size = 12))


# load dataset
data <- read_csv("HLWconjoint20210330.csv")


# processing of data
data2 <- select(data, no, time, rl, group, hlw_su, hlw_fsu, 
                hlw_dis, hlw_pol, hlw_so, hlw_ad, hlw_pop, hlw_inc)
data2$hlw_fsu <- ifelse(data2$hlw_fsu == data2$rl, 1, 0)
data2 <- mutate(data2, 
                `Distance` = recode(hlw_dis, `1` = "less than 1 km", 
                                    `2` = "1-5 km", `3` = "5-20 km", 
                                    `4` = "20 km or more"), 
                `Politics` = recode(hlw_pol, `1` = "job creation", 
                                    `2` = "public services", 
                                    `3` = "local economy", 
                                    `4` = "society"), 
                `Social` = recode(hlw_so, `1` = "70% agreed, 30% disagreed", 
                                  `2` = "50% agreed, 50% disagreed", 
                                  `3` = "30% agreed, 70% disagreed"), 
                `Administration` = recode(hlw_ad, `1` = "none", 
                                          `2` = "childbirth and childcare", 
                                          `3` = "elderly care", `4` = "education", 
                                          `5` = "public facility", `6` = "water charge"), 
                `Population` = recode(hlw_pop, `1` = "19,000", 
                                      `2` = "20,000", `3` = "21,000"), 
                `Income` = recode(hlw_inc, `1` = "3.8 million yen", 
                                  `2` = "4 million yen", `3` = "4.2 million yen"))
data2$Distance <- factor(
  data2$Distance, 
  levels = c("less than 1 km", "1-5 km", "5-20 km", "20 km or more")
)
data2$Politics <- factor(
  data2$Politics, 
  levels = c("society", "job creation", "public services", "local economy")
)
data2$Social <- factor(
  data2$Social, 
  levels = c("30% agreed, 70% disagreed", "50% agreed, 50% disagreed", "70% agreed, 30% disagreed")
)
data2$Administration <- factor(
  data2$Administration, 
  levels = c("none", "childbirth and childcare", "elderly care", 
             "education", "public facility", "water charge")
)
data2$Population <- factor(
  data2$Population, 
  levels = c("19,000", "20,000", "21,000")
)
data2$Income <- factor(
  data2$Income, 
  levels = c("3.8 million yen", "4 million yen", "4.2 million yen")
)


# AMCE - forced choice
res1 <- cj(
  data2, 
  hlw_fsu ~ Distance + Politics + Social + Administration + Population + Income, 
  id = ~ no, estimate = "amce", level_order = "descending"
)
res1$level <- as.character(res1$level)
res1 <- rbind(
  res1, 
  c("hlw_fsu", "amce", "Distance", "(Distance attribute)", NA, NA, NA, NA, NA, NA), 
  c("hlw_fsu", "amce", "Politics", "(Political attribute)", NA, NA, NA, NA, NA, NA), 
  c("hlw_fsu", "amce", "Social", "(Social attribute)", NA, NA, NA, NA, NA, NA), 
  c("hlw_fsu", "amce", "Administration", "(Administrative attribute)", NA, NA, NA, NA, NA, NA), 
  c("hlw_fsu", "amce", "Population", "(Population attribute)", NA, NA, NA, NA, NA, NA), 
  c("hlw_fsu", "amce", "Income", "(Average income attribute)", NA, NA, NA, NA, NA, NA)
)
res1$level <- factor(
  res1$level, 
  levels = c(
    "20 km or more", "5-20 km", "1-5 km", "less than 1 km", "(Distance attribute)", 
    "local economy", "public services", "job creation", "society", "(Political attribute)", 
    "70% agreed, 30% disagreed", "50% agreed, 50% disagreed", "30% agreed, 70% disagreed", "(Social attribute)", 
    "water charge", "public facility", "education", "elderly care", "childbirth and childcare", "none", "(Administrative attribute)",  
    "21,000", "20,000", "19,000", "(Population attribute)", 
    "4.2 million yen", "4 million yen", "3.8 million yen", "(Average income attribute)"
  )
)
res1$estimate <- as.numeric(res1$estimate)
res1$std.error <- as.numeric(res1$std.error)
res1$z <- as.numeric(res1$z)
res1$p <- as.numeric(res1$p)
res1$lower <- as.numeric(res1$lower)
res1$upper <- as.numeric(res1$upper)
res1.plot <- ggplot(
  res1, aes(x = reorder(level, as.numeric(feature) * -1), y = estimate, 
            ymax = upper, ymin = lower, shape = feature)) + 
  geom_hline(yintercept = 0, color = "gray60") + 
  geom_pointrange(size = 0.2) + coord_flip() + ylim(-0.35, 0.35) + 
  labs(y = "Estimated AMCE", title = "Forced Choice") + 
  theme_minimal(base_size = 12) + 
  theme(plot.title = element_text(hjust = 0.5), 
        axis.title.y = element_blank(), 
        legend.position = "none", 
        legend.title = element_blank())
plot(res1.plot)


# AMCE - rating
res2 <- cj(
  data2, 
  hlw_su ~ Distance + Politics + Social + Administration + Population + Income, 
  id = ~ no, estimate = "amce", level_order = "descending"
)
res2$level <- as.character(res2$level)
res2 <- rbind(
  res2, 
  c("hlw_su", "amce", "Distance", "(Distance attribute)", NA, NA, NA, NA, NA, NA), 
  c("hlw_su", "amce", "Politics", "(Political attribute)", NA, NA, NA, NA, NA, NA), 
  c("hlw_su", "amce", "Social", "(Social attribute)", NA, NA, NA, NA, NA, NA), 
  c("hlw_su", "amce", "Administration", "(Administrative attribute)", NA, NA, NA, NA, NA, NA), 
  c("hlw_su", "amce", "Population", "(Population attribute)", NA, NA, NA, NA, NA, NA), 
  c("hlw_su", "amce", "Income", "(Average income attribute)", NA, NA, NA, NA, NA, NA)
)
res2$level <- factor(
  res2$level, 
  levels = c(
    "20 km or more", "5-20 km", "1-5 km", "less than 1 km", "(Distance attribute)", 
    "local economy", "public services", "job creation", "society", "(Political attribute)", 
    "70% agreed, 30% disagreed", "50% agreed, 50% disagreed", "30% agreed, 70% disagreed", "(Social attribute)", 
    "water charge", "public facility", "education", "elderly care", "childbirth and childcare", "none", "(Administrative attribute)",  
    "21,000", "20,000", "19,000", "(Population attribute)", 
    "4.2 million yen", "4 million yen", "3.8 million yen", "(Average income attribute)"
  )
)
res2$estimate <- as.numeric(res2$estimate)
res2$std.error <- as.numeric(res2$std.error)
res2$z <- as.numeric(res2$z)
res2$p <- as.numeric(res2$p)
res2$lower <- as.numeric(res2$lower)
res2$upper <- as.numeric(res2$upper)
res2.plot <- ggplot(
  res2, aes(x = reorder(level, as.numeric(feature) * -1), y = estimate, 
            ymax = upper, ymin = lower, shape = feature)) + 
  geom_hline(yintercept = 0, color = "gray60") + 
  geom_pointrange(size = 0.2) + coord_flip() + ylim(-0.4, 0.4) + 
  labs(y = "Estimated AMCE", title = "Rating") + 
  theme_minimal(base_size = 12) + 
  theme(plot.title = element_text(hjust = 0.5), 
        axis.title.y = element_blank(), 
        legend.position = "none", 
        legend.title = element_blank())
plot(res2.plot)


# MM - forced choice
res3 <- cj(
  data2, 
  hlw_fsu ~ Distance + Politics + Social + Administration + Population + Income, 
  id = ~ no, estimate = "mm", level_order = "descending"
)
res3$level <- as.character(res3$level)
res3 <- rbind(
  res3, 
  c("hlw_fsu", "mm", "Distance", "(Distance attribute)", NA, NA, NA, NA, NA, NA), 
  c("hlw_fsu", "mm", "Politics", "(Political attribute)", NA, NA, NA, NA, NA, NA), 
  c("hlw_fsu", "mm", "Social", "(Social attribute)", NA, NA, NA, NA, NA, NA), 
  c("hlw_fsu", "mm", "Administration", "(Administrative attribute)", NA, NA, NA, NA, NA, NA), 
  c("hlw_fsu", "mm", "Population", "(Population attribute)", NA, NA, NA, NA, NA, NA), 
  c("hlw_fsu", "mm", "Income", "(Average income attribute)", NA, NA, NA, NA, NA, NA)
)
res3$level <- factor(
  res3$level, 
  levels = c(
    "20 km or more", "5-20 km", "1-5 km", "less than 1 km", "(Distance attribute)", 
    "local economy", "public services", "job creation", "society", "(Political attribute)", 
    "70% agreed, 30% disagreed", "50% agreed, 50% disagreed", "30% agreed, 70% disagreed", "(Social attribute)", 
    "water charge", "public facility", "education", "elderly care", "childbirth and childcare", "none", "(Administrative attribute)",  
    "21,000", "20,000", "19,000", "(Population attribute)", 
    "4.2 million yen", "4 million yen", "3.8 million yen", "(Average income attribute)"
  )
)
res3$estimate <- as.numeric(res3$estimate)
res3$std.error <- as.numeric(res3$std.error)
res3$z <- as.numeric(res3$z)
res3$p <- as.numeric(res3$p)
res3$lower <- as.numeric(res3$lower)
res3$upper <- as.numeric(res3$upper)
res3.plot <- ggplot(
  res3, aes(x = reorder(level, as.numeric(feature) * -1), y = estimate, 
            ymax = upper, ymin = lower, shape = feature)) + 
  geom_hline(yintercept = 0.5, color = "gray60") + 
  geom_pointrange(size = 0.2) + coord_flip() + ylim(0.25, 0.75) + 
  labs(y = "Estimated Marginal Mean", title = "Forced Choice") + 
  theme_minimal(base_size = 12) + 
  theme(plot.title = element_text(hjust = 0.5), 
        axis.title.y = element_blank(), 
        legend.position = "none", 
        legend.title = element_blank())
plot(res3.plot)


# MM - rating
res4 <- cj(
  data2, 
  hlw_su ~ Distance + Politics + Social + Administration + Population + Income, 
  id = ~ no, estimate = "mm", level_order = "descending"
)
res4$level <- as.character(res4$level)
res4 <- rbind(
  res4, 
  c("hlw_su", "mm", "Distance", "(Distance attribute)", NA, NA, NA, NA, NA, NA), 
  c("hlw_su", "mm", "Politics", "(Political attribute)", NA, NA, NA, NA, NA, NA), 
  c("hlw_su", "mm", "Social", "(Social attribute)", NA, NA, NA, NA, NA, NA), 
  c("hlw_su", "mm", "Administration", "(Administrative attribute)", NA, NA, NA, NA, NA, NA), 
  c("hlw_su", "mm", "Population", "(Population attribute)", NA, NA, NA, NA, NA, NA), 
  c("hlw_su", "mm", "Income", "(Average income attribute)", NA, NA, NA, NA, NA, NA)
)
res4$level <- factor(
  res4$level, 
  levels = c(
    "20 km or more", "5-20 km", "1-5 km", "less than 1 km", "(Distance attribute)", 
    "local economy", "public services", "job creation", "society", "(Political attribute)", 
    "70% agreed, 30% disagreed", "50% agreed, 50% disagreed", "30% agreed, 70% disagreed", "(Social attribute)", 
    "water charge", "public facility", "education", "elderly care", "childbirth and childcare", "none", "(Administrative attribute)",  
    "21,000", "20,000", "19,000", "(Population attribute)", 
    "4.2 million yen", "4 million yen", "3.8 million yen", "(Average income attribute)"
  )
)
res4$estimate <- as.numeric(res4$estimate)
res4$std.error <- as.numeric(res4$std.error)
res4$z <- as.numeric(res4$z)
res4$p <- as.numeric(res4$p)
res4$lower <- as.numeric(res4$lower)
res4$upper <- as.numeric(res4$upper)
res4.plot <- ggplot(
  res4, aes(x = reorder(level, as.numeric(feature) * -1), y = estimate, 
            ymax = upper, ymin = lower, shape = feature)) + 
  geom_pointrange(size = 0.2) + coord_flip() + ylim(2.5, 3.5) + 
  labs(y = "Estimated Marginal Mean", title = "Rating") + 
  theme_minimal(base_size = 12) + 
  theme(plot.title = element_text(hjust = 0.5), 
        axis.title.y = element_blank(), 
        legend.position = "none", 
        legend.title = element_blank())
plot(res4.plot)
